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ABSTRACT 

We present a new exact method to numerically compute the thermodynamical properties of 
an interacting Bose gas in the canonical ensemble. As in our previous paper (Phys. Rev. 
A 63 023606 (2001) ), we write the density operator p as an average of Hartree dyadics 
\N : (j)i){N : and we find stochastic evolution equations for the wave functions 0i^2 
such that the exact imaginary-time evolution of p is recovered after average over noise. In 
this way, the thermal equilibrium density operator can be obtained for any temperature 
T. The method is then applied to study the thermodynamical properties of a homogeneous 
one-dimensional A-boson system: although Bose-Einstein condensation can not occur in the 
thermodynamical limit, a macroscopic occupation of the lowest mode of a finite system is 
observed at sufficiently low temperatures. If ksT 3> p, the main effect of interactions is to 
suppress density fluctuations and to reduce their correlation length. Different effects such 
as a spatial antibunching of the atoms are predicted for the opposite fc^T < p regime. Our 
exact stochastic calculations have been compared to existing approximate theories. 

PACS: 03.75.Fi, 05.10. Gg. 42.50.-p 



I. INTRODUCTION 

The recent advances in the production and manipulation of ultra-cold atomic samples have 
opened the way to the investigation of the thermodynamical and dynamical properties of 
interacting Bose gases in different geometries and dimensionalities |]T]-§[; in three dimensional 
traps, Bose-Einstein condensation was realized a few years ago |^ and now there is a great 
deal of experimental interest in lower- dimensional systems, for which a remarkably different 
physics is predicted [p|-p!0|]. 

In presence of a condensate, i.e. when the coherence of the Bose field extends throughout 
the whole sample, a mean-field theory can be applied, which completely neglects the non- 
condensed atoms and describes the system in terms of the classical field corresponding to the 
condensate wave function O]. In three dimensions, such an approach is clearly valid only 



when a small non-condensed fraction is present, which implies a low temperature and weak 
interaction regime; for low- dimensional systems a further constraint has to be imposed on the 
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spatial size of the system, which has to be smaller than the finite coherence length of the field. 
Within mean-field theory, weak thermal excitations on top of a well populated condensate 
can be taken into account by means of the Bogoliubov approach which describes the low-lying 
part of the many-body spectrum as a collection of uncoupled harmonic oscillators |1^-[15|. 
Clearly, this approximated theory is valid only at relatively low temperatures, when the 
density of thermal excitations is small enough for their mutual interaction to be neglected. 

Exact results for the thermodynamical properties of interacting Bose gases at any tempera- 
ture can in principle be obtained using Quantum Monte Carlo techniques [16|, but such an 
approach is subjected to strict limitations in its applicability domain: first of all, it requires 
that the matrix elements of the Hamiltonian are real in the position representation; this 
condition excludes rotating systems such as the ones used to study vortices. The privileged 
role attributed to the position representation makes some observables hard to be sampled, 
e.g. the momentum distribution function or the occupation number statistics of a given 
mode. Finally, the thermodynamical calculations performed with the quantum Monte Carlo 
method can not be prolongated to real time so to determine the dynamical properties of the 
system. 

In the present paper we present a new exact and numerically tractable approach to the 
calculation of thermodynamical properties of an interacting Bose gas in the canonical en- 
semble. This computational scheme is not restricted to Hamiltonians with real matrix 
elements in real space and gives access to all observables of the system. Also, this computa- 



tional scheme can be easily integrated with our reformulation of the dynamical evolution ||T7 
and therefore used to describe the response of finite-temperature samples to perturbations 
away from equilibrium. The method is based on the stochastic evolution of Hartree dyadics 
(T = |iV : 01 )(iV : 02| where 0i_2 are non- normalized wave functions. In our previous pa- 
per [0 we showed that it is possible to choose the stochastic evolution of 0i,2 in such a 
way that the average evolution of o over the stochastic noise exactly reproduces the full 
many-body time-evolution. Here, we show that similar arguments can be used to find a 
stochastic evolution which exactly reproduces the imaginary-time evolution of the Bose gas 
and therefore can be used to determine the thermodynamical properties of the interact- 
ing A^-boson system in the canonical ensemble. Differently from previous theories such as 
Positive-P representation, our scheme is not subjected to instability problems [|18 



After the general discussion on the method, we shall report its first application to a homoge- 
neous one-dimensional interacting Bose gas: such a model is expected to accurately describe 
the physics of Bose condensates tightly confined in the transverse dimensions so to create 
an effective one-dimensional geometry 

The outline of the paper is the following: in Sec.|| we present the general theory underlying 
the method and we discuss the stochastic equations for the Hartree wave functions 0i,2 to 
be used for the study of dynamics and thermodynamics of the interacting Bose gas. 

In Sec.|T|, we give details on the algorithm used in numerical simulations in order to enhance 
the efficiency of Monte Carlo calculations; in particular, a large reduction in the statistical 
noise can be obtained by means of an importance sampling with an appropriate a •priori 
distribution function. 

In Sec.[V|, we shall apply the method to the thermodynamical properties of a homogeneous 
one-dimensional interacting Bose gas in the canonical ensemble and we shall compare the 
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results of our stochastic calculations with existing approximate theories. 
Conclusions are finally given in Sec.0 



II. THE STOCHASTIC WAVE FUNCTION APPROACH 



The Hamiltonian of a trapped interacting Bose gas in D dimensions can be written in a 
second-quantized form as 



H = y"d^x^t(x)/io^(x) + l J y"(i^x(iV^t(x)^t(x')l^„,(x-x')^(x')^(x) 



:i) 



where Hq = ~|^V^ + Kxtlx) is the single particle Hamiltonian in the external confining 
potential \4xt(x). Interactions are modeled by the two-body potential V^nt(x — x') which 
is assumed to be repulsive. The field operators ^I/(x) satisfy Bose commutation relations 
[^(x),^t(x')] = 5(x-x'). 

Although the exact dynamics of the many-body system could in principle be obtained from 
the quantum mechanical equation of motion for the iV-body density operator p 

any explicit calculation results impossible for the experimentally relevant multi-mode sys- 
tems because of the huge dimensionality of the corresponding A^-body Hilbert space. 

In this section we shall discuss a reformulation of the bosonic many-body problem in terms 
of stochastic wave functions. As we have shown in the A^-body density operator of the 



system can always be written as a statistical average of Hartree operators \N : (j)i){N : (j)2\ 
in which all N atoms share the same, not necessarily normalized wave function (Sec . [II A]) . In 
the same paper, we have shown that the exact time-evolution of the full many-body system 
can be reproduced by an appropriately chosen stochastic evolution of the pair of single- 
particle wave functions 0i,2 (Sec. |ll B| ). The central theoretical result of the present paper 
is presented in Sec. |ll Q a similar stochastic evolution of the single particle wave functions 
01 2 is able to reproduce the imaginary-time evolution: 

^ = m, P(r)} = -\ {n p{r) + p(r) H) (3) 

of the many body system with initial condition p(r = 0) = Id; this gives access to the 
thermodynamical properties of the interacting Bose gas in the canonical ensemble at any 
desired temperature T as the solution of (||) after an evolution "time" r = l/Zc^T is p(r = 



A. The Fock state Hartree dyadic ansatz 



As in our precedent paper |]T^, we consider the A^-body Fock state Hartree ansatz 

a= |Ar:0i)(iV:02| (4) 
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in which the two Hartree wave functions 0i,2 in the bra and the ket are in general different 
and not normalized. 

Given the remarkable identity denotes the functional integration over the set of 

wave functions of unit norm \\(f)\\ = 1) 

ld = jv(t)\N : (t)){N : (t)l (5) 

any A^-body density operator p can be expanded over the family of dyadics (H) 

p = J y"p0iD02P(0i,02)|iV:0i)(iV:02| (6) 

with a positive probability distribution P(0i,02)- Note that this expansion is by no means 
unique and several distinct V{(f)i,4>2) can be found which correspond to the same physical 
density operator p. 

From the expansion , the mean values of the observables can be written as mean values of 
combinations of 0i and 02 over the probability distribution P(0i,02)- As simple examples, 
the one-body density matrix p^^'^ (x, x') can be written as 

p«(x,x') = (>I/t(x')^(x)) = iV0i(x)0^(xO(02|0i)^-^ (7) 

and the second-order correlation function C^^''(x, x') as 

C(2)(x,x') = {^^^(x)^^(x')^(x')^(x)^ = 

= N{N - 1) 0i(x)0i(xO0^(x')0^(x)(02|0i)^-2. (8) 



B. Real-time evolution: dynamics 



The central result of the paper [|T3 was the proof of the possibility of reproducing the exact 
time evolution of p by means of a Ito stochastic evolution for 0i 2 

(j)a{t + dt) = (j)ait) + Fa dt + dB^ (a = 1,2); (9) 



in all the paper, the noise term dBa is treated in the standard Ito formalism [|I^: it has 
a zero mean dB^ and dB"^ oc dt; the deterministic contribution is given by the force term 
Fadt. 

In the infinitesimal time interval [t,t + dt], the two wave functions 0i(x) and 02 (x) are 
assumed to evolve according to Ito stochastic differential equations (P). After dt, the dyadic 
a{t) = : 01 )( iV : 02| has therefore evolved into 

a{t + dt) = \N : 0i + rf0i)(iV : 02 + rf02| (10) 

where c?0i^2 are defined according to (y). Imposing that this evolution reproduces in the 
average the exact many-body evolution 

a{t + dt) = a{t) + '^{'Hait) - a{t)n) (11) 
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gives a set of constraints on the explicit form of the stochastic differential equations (H). 

As we have shown in many stochastic evolutions can be found which all satisfy the 
exactness constraints: their stability and statistical properties are however significantly dif- 
ferent. In the following we shall limit our attention to the so-called simple scheme minimizing 
the growth of the statistical variance of a: thanks to its stability, this scheme has in fact 
been found to be the most efficient for practical simulations. 

Explicitly, the stochastic evolution of the wave functions 0i 2 reads 



ih 



+ K.M) + / rfx' l^„t(x - x') |0«(x')|' + 

2"^ ||0a|| J 

(iV- 1) (0<,0«|Mnt|0a0a) 



, dt + dB^] (12) 



as expected, the deterministic term closely resembles the Gross-Pitaevskii evolution of mean- 
field theory; the noises dBi and dB2 are statistically independent and are determined by the 
correlation properties 



(i5«(x) dB^{^) = %QlQt [\/int(x - X')0a(x)04x')] 

in 



[a 



1,2) 



(13) 



where is the projector onto the subspace orthogonal to 0q(x). Thanks to the fact 
that the deterministic term is norm-preserving, a mathematically well-defined solution ex- 
ist for all times so that the simulation does not suffer from the divergence prob- 
lems encountered by previous stochastic approaches such as Gardiner's Positive-P distri- 
bution [|T^ , pO| -|22[| . The statistical error, although exponentially growing with time propor- 
tionally to exp[2A^V^nt(0)t/^], is always finite Q : this guarantees that rehable physical results 
can be obtained by means of a Monte Carlo simulation of the stochastic differential equa- 
tions; the norm-preserving character of the deterministic force ensures that no intrinsically 



spiking trajectories ||T8[ can occur in the numerical calculations . This fundamental fact has 
been verified in Monte Carlo simulations of simple systems. 

Since any A^-body density operator p{t) can be expanded in Hartree dyadics (^, the time 
evolution of the many body system can be reproduced by the stochastic approach starting 
from any initial condition. In |l7| we used as initial state a fully coherent mean-field state 
p = \N : (j)) [N : (j)\\'m order to be able to get information on the full dynamics keeping track 
of all quantum effects, a more physical initial state has to be chosen. In the following part of 
the paper we shall show how a related stochastic wave function approach can be found which 
is able to sample the thermal equilibrium density operator peq{(3) for a many-body system 
of N interacting bosons at a given (3 = l/ksT in the canonical ensemble. A simulation of 
the dynamics of the many body system at finite temperature f3 can be then performed on 
the base of (H) using PeqiP) as the initial state. 



^ This result holds in the case of repulsive interaction potentials with a positive Fourier trans- 

int(x) 



form 0, for which yint(O) > |M 
vanishes in the absence of interactions 



This means, in particular, that the statistical error only 
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C. Imaginary-time evolution: thermodynamics 



It is a well-known fact of quantum statistical mechanics that the unnormalized thermal 
equilibrium density operator at a given temperature T in the canonical ensemble can be 
written as 

Pe,{(3) = e-^^ (14) 
with (3 = 1/kBT. Starting from the infinite temperature equilibrium state 

p,,(/3 = 0)=Id, (15) 
the density operator at any given jS can be obtained by means of the imaginary-time evolution 

dpeqir) 



dr 



{H, Pe,(r)} = {Hpegir) + Pe,{T) H) 



(16) 



during a "time" (3. The mean values of the observables are then easily obtained from the 
equilibrium density operator 



1 



Tr [p 



fTr [pe,0]. 



(17) 



Given the similarity of the real- and imaginary-time evolution equations (0) and (|TB|), the 
same arguments which led to (|1^,0) can be generalized to the imaginary-time evolution: 
the stochastic equations for the wave functions 0i^2 within the simple scheme now read 



dr 
T 



+ K..(x) + ^^^—^ / dx' VU^ - x'l 



2m 



(AT - 1) (0„0,|Vlnt|0a0a) 



:i8) 



and the correlation functions of the statistically independent noise terms dBi 2 are fixed by 

dr 



dB^i^) dB4^') = -^QlQa [^nt(x - x')0a(x)0«(x')] ; 



(19) 



it is apparent how these equations are obtained from the real-time ones by means of the 
simple substitution dt —ih^. 

As the starting point of the imaginary-time evolution, one has to choose a representation of 
the identity operator describing the T = oo equilibrium state; even if not the most efficient 
one, the simplest choice is clearly the one dictated by (1^) 



Pe,(r = 0)= V4>\N:<P){N:(P\; 



(20) 



in the following Sec.|T|, we shall discuss how the efficiency of an actual Monte Carlo simu- 
lation can be sensibly improved by choosing a non-trivial a priori probability distribution. 
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Differently from what happened for the real-time evolution, the deterministic term now 
resembles an imaginary-time Gross-Pitaevskii equation and does not preserve any more the 
norm ||0i,2|| of the wave functions. In particular, the deterministic contribution to the norm 
variation is equal to 

dUX = -eGp[(t>a\Ua\? dr (21) 
where the mean-field energy eGp[0] is defined in terms of the normalized wave function 

(t^a = 4>a/\\4>a\\ aS 

(^Gp[4>a] = (0a I ho \(f)a) + ^ V^nt \ (pa4>a) ■ (22) 

Since interactions are assumed to be repulsive and the atoms are assumed to be trapped 
(the spectrum of Hq is bounded from below) the mean- field energy ecp is also bounded 
from below; this guarantees that the stochastic equations have regular and non-exploding 
solutions defined at all r's and the stochastic dynamics can be safely simulated by means of 
Monte Carlo techniques. 



III. THE SIMULATIONS 



The mathematical manipulations of the previous section have led to a reformulation of the 
interacting A^-body problem in terms of a pair of simple stochastic differential equations for 
the wave functions 0i_2- This is a good starting point for a Monte Carlo simulation of the full 
quantum dynamics of the many body system: in we presented first examples of real-time 
simulations, in the present paper we shall present imaginary-time ones based on the theory 
of Sec.|ll C|. The present section is devoted to a general analysis of the simulation algorithm. 



while the next Sec.|3 will describe the application of the method to the thermodynamical 
properties of a one- dimensional interacting Bose gas. 

The simplest Monte Carlo implementation of the simulation scheme consists in randomly 
choosing an initial wave function (p on the unit sphere with an uniform distribution as done in 
(PPP; the two wave functions 0i^2 are then evolved from their initial value (j)i^2 = 'P according 
to the stochastic evolution (|18D including the noise term (pj]). 

Unfortunately, this simple procedure is not efficient to calculate the expectation value of the 
observables at a given (3. Since the largest fraction of the randomly chosen wave functions 
initially have a large mean- field energy (0), their deterministic evolution according to (^I]) 
gives a fast reduction of the norm. The explicit expression of the physically interesting 
observables, (|^) or (^, contains the norm of the wave function raised to a large power of 
the order A^; in particular, the trace of the denominator of (|I^) involves a sum of stochastic 
realizations of 

(02i0i)^cxii0iini02ir. (23) 

The contribution of most realizations to this sum is therefore negligible so that a large 
fraction of the computational time is effectively wasted and the convergence of the stochastic 
average is impractically slow. 
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A technique which is currently used in Monte Carlo integration to overcome this kind of 
problems is the so-called importance sampling, a general and comprehensive discussion of 
which can be found in §7.8]. In our specific case, such a technique is equivalent to using 
the following equivalent representation of the identity operator instead of the simple (|^) 

P.,(. ^ 0) ^ jmv, (24) 

the initial wave function is now chosen according to the (arbitrarily chosen) a priori proba- 
bility distribution P[4>], while the weight of the realization is correspondingly divided by a 
factor P[(f)], which amounts to a global rescaling of 0. By choosing a P[0] concentrated on 
the realizations with largest weight at "time" f3, an effectively flatter integrand is obtained 
for the expectation values of the observables which means a more efficient Monte Carlo sam- 
pling. Although different observables lead to different optimization schemes, a reasonable 
choice if we are mainly interested in few-body observables such as spatial and momentum 
densities and low-order correlation functions can be an optimization scheme taking the trace 
of p as an optimization criterion: we wish that each dyadic : 4'i){N : 02| in the simulation 
has approximately the same trace (|2|). 

For a non-interacting gas, no stochastic noise is present in the evolution equation (|TBp, so 
the two wave functions are always equal 0i = 02- Let us first consider a dyadic a with initial 
value : 0(0) )( : 0(O)|/P[0(O)], where 0(0) is distributed over the unit sphere with the 
probability density P[0(O)]. After an imaginary-time evolution during r, the trace of a is 
given by ||0(r)|p^/P[0(O)] where 0(r) is the result of the imaginary-time evolution starting 
from 0(0). The imaginary-time evolution in the absence of interactions is simply written in 
the basis of eigenmodes of the trapping potential, that is the eigenvectors of h^. We label 
these eigenmodes with k and we note their eigenenergy Ek: 

Mr) = 0.(O)e--^^/2; (25) 

0fc(r) denotes the component of the wave function 0(r) on the eigenmode k. In this simple 
case, a unit trace for the dyadic a at "time" t = P can be obtained for the following choice 
of Pf 




^[0(0)] = wmr = [J2\^km'e-^^>' . (26) 



Even if not the optimal one, this choice is able to provide good efficiency also in the inter- 
acting case as we shall see on numerical examples. 

In the next subsections we shall discuss two different methods to numerically generate ran- 
dom initial wave functions according to the probability distribution P[0] of ([26|) . The first 
method (sec. [Ill A]) has to be used when no condensate is expected to be present. The second 



one (sec. |IIIB ) is to be preferred at low temperatures when there is a macroscopic occupation 
of a single mode; in this range of parameters, in fact, this second method is able to provide 
a quicker and more accurate sampling of P[0], but it is extremely slow in the absence of a 
condensate. 
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A. Brownian motion sampling 



A simple way to generate a random variable with a given probability distribution function 
P[0] consists in constructing a Brownian motion with a deterministic drift term F^'^'^*) [0] 
and a Langevin noise term dcj)^^*'^: 

d4> = F^'^^'\(f)]dt + d4>^''^ (27) 

whose equilibrium distribution function is exactly P[(f)]; if the random variable is evolved 
for a sufficiently long time tmax, its value will be independent from the initial value and 
distributed according to -P[0]. In practical simulations, tmax is chosen in a heuristic way 
so to make results not to change if it is further increased. Note that t is here a fictitious 
evolution time that we take dimensionless. 

Provided the force term is integrable 

= -2V^.f/[0] (28) 

and the diffusion is constant and isotropic 

d(t>i''^d<l)l\''^ =26k,k'dt, (29) 
d</,f ) = (30) 

the equilibrium density of the Brownian motion (|27|) can be shown to be equal to 

P,J0] = e-2^[^l. (31) 

In our specific case, the Brownian motion has to be confined on the ||0|| = 1 sphere: this 
is actually done by keeping only the projection of both the deterministic force and the 
stochastic noise along the sphere surface (i.e. orthogonal to 0) and then renormalizing the 
wave function to unity after each time step; the renormalization procedure is equivalent to 
adding a deterministic force along (p which keeps the wave function on the sphere. 

Explicitly, the deterministic and stochastic increments can be written in terms of the the 
projector = 1 - as 



d<P'st = Qi (32) 
FLl^] = Qi F^'''\<f>] (33) 

with the force F^*^*^*) being given by the derivative with respect to 0^(0) of the logarithm of 

- V- 1^ |2e-/3s,,- (34) 



In order to avoid systematic errors, an actual calculation requires a temporal step dt such that 
both the deterministic F^g^ dt and the stochastic d(j)'g^ variations of are small as compared 
to II 011 = 1, which respectively means N dt 1 and N" dt <^ 1, Af being the number of 
modes actually used in the simulation. 
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B. Bogoliubov rejection sampling 



While the samphng method discussed in the previous section can be used for any choice 
of physical parameters, in the presence of a large population in a single mode, a different 
method based on a Bogoliubov sampling can be more efficient; it is also free from arbitrary 
parameters like tmax or dt of the previous Brownian motion sampling which are sources of 
systematic errors. 

Taking into account the wave function normalization ||0|| = 1, we can eliminate from the 
probability distribution (^) the field variable 0o using |0oP = 1 — J2kf^o l^fcP (its phase is 
uniformly distributed) and derive a probability density for the 0a;'s, A; 7^ 0, in terms of the 
(j)kjLo only: 



TV 



dAf 



TC 



E 






















e 


1 




^k J 




k^O 



(35) 

ky^Q 



where the d(f)kS stand for the double integral over the real and the imaginary parts of 0^; 
the energy Eq of the ground k = mode is set to Eq = 0. The coefficients Afc^o are defined 
as 



Afe = (1 



(36) 



and satisfy > 1. Up to a global multiplicative factor, the distribution function (35) can 
be rewritten in terms of the reduced field ipk = <Pk/ v^A^ as 



dAf- 




1 - E Afcl^/; 



kj^O 



(37) 



and in this form it is easily sampled: ipk^o are generated with a spherically symmetric 
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probability distribution proportional to 



(40) 

and then all the realizations which do not satisfy the G function in ( pT]) are rejected i.e. 
the ones for which 

$^Afc|V'fcP>l. (41) 

At the end of this procedure, (p^s are obtained which are distributed according to the desired 
distribution law 

The efficiency of this method clearly depends on the fraction of realizations which have to 
be rejected at some stage; the step which strongly depends on the physical parameters is 
the rejection governed by condition (|T): the condition (|T) can be reformulated in terms 



of the number = N\(f)k\ of particles in the k excited mode as 

J2Nk>N; (42) 

the larger the occupation of the ground A; = mode, the smaller the number of realizations 
rejected at this stage and thus the quicker the sampling. In presence of a macroscopic 
occupation, the most probable occupation number of the k = mode is non-zero (cfr. 
Sec. [lV B"3D so that the rejection occurs very rarely and the Bogoliubov sampling method is 



very rapid. 



IV. APPLICATION: THERMODYNAMICAL PROPERTIES OF A ID 

HOMOGENEOUS BOSE GAS 

In the previous sections we have described the principles of the stochastic approach to 
the imaginary-time evolution of an interacting Bose gas as well as the crucial points of its 



^Thanks to the inequality 

e-^E.^ol'^.P>(l_^|^,|2)^. (38) 

the distribution function (^) can be sampled by means of the rejection niettiod §7.3] with the 
Gaussian as the comparison function: this latter is sampled by standard methods p3| , §7.2]; for 
each reahzation, we then choose a random number z uniformly distributed in [0, 1]: if 



\2\N 



^> _w^,^i,/„j2 (39) 



the reahzation is rejected, otherwise it is retained. The remaining reahzations are automatically 
distributed according to (^). 
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actual implementation in a Monte Carlo simulation. The present section is devoted to the 
application of the method developed in the previous sections to a physically relevant system 
such as interacting bosons in a one dimensional box of finite size L with periodic boundary 
conditions. 

The one-dimensional regime is not far from reach in present experiments on transversally 
confined Bose condensates In particular, the phase fluctuation of the field when the 
coherence length is smaller that the spatial extension of the gas have been recently observed 
in a harmonic potential . 



A. Ideal Bose gas in a large box: quantum degeneracy 

It is a well-known fact of thermodynamics that Bose-Einstein condensation does not occur 
in a one-dimensional and non-interacting homogeneous system in the thermodynamic limit 
L — > oo at a constant density n = N/L |jl3|. Even in the absence of Bose condensation, 
effects coming from the Bose statistics are however apparent in the quantum degenerate 
regime 



24 



T < Tdeg = 43 

TUKb 



i.e. when the thermal de Broglie wavelength 



is larger than the mean inter particle spacing d = n~^; in this case, Bose statistics strongly 
enhances population of the lowest modes so that the characteristic field coherence length 

4=4^ = ^ (45) 
^ mkeT 27r ^ ^ 

turns out to be much longer than in the case of a non-degenerate system for which the 
coherence length is equal to \th/V^- 

In fig.|I] the results of our stochastic simulations for a non-interacting gas at temperatures 
both above and below the degeneracy temperature are compared with the analytical predic- 
tions of the grand-canonical ensemble. Provided the size of the box L is chosen sufficiently 
larger than the coherence length ic (which is the case of fig.|I|); excellent agreement is found; 
in this limit the canonical and the grand-canonical ensembles are in fact equivalent. In the 
next section specific effects due to a fixed particle number and the finite size of the system 
will be addressed. 

Differently from the 3D case in which a true condensate is present in the thermodynamical 
limit, the ID first-order correlation function (fig.|^) 

gW{x) = -{^\x)^{0)) (46) 

Th 
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has a vanishing long distance hmit g^^\x oo) = and a similar behavior is found in the 
second-order correlation function (fig.^): 

gi^)(x) = i-(^U0WUx)^(x)^(0)y, (47) 

in the short distance limit, g^"^^ tends to the value 2 typical of thermal non-interacting Bose 
systems; for large x, g^'^\x) exponentially decays to the value 1 with a decay constant equal 
to i'c/2: the characteristic correlation length scales for the phase ( ^61) and density ( ^Tf ) 
fluctuations are therefore of the same order. 

All the calculations presented in this section have been performed by means of the Brownian 
motion sampling discussed in Sec. [lII A] ; in the present L 3> £c regime, the efficiency of a 
Bogoliubov sampling would be definitely poorer. 



B. Ideal Bose gas in a small box: Bose-Einstein condensation 

In the present section we shall discuss the case of a true Bose-Einstein condensate in a finite 
ID box: in particular, we shall study the effect of a macroscopic occupation of a single mode 
on different observables; differently from other techniques , our method is in fact able to 
give access to all the observables of the system. 



1. Macroscopic occupation of a single mode 

In the previous section, we have shown that the phase coherence in a ID sample extends 
only over a finite length ic and the absence of a true condensate in the thermodynamical 
limit corresponds to a vanishing limit g^^\x oo). On the other hand, if the sample we 
are considering has a length L comparable to or smaller than the coherence length ic, phase 
coherence extends over the whole sample. 

According to the Wiener- Khint chine theorem [PT|, the momentum distribution function A^(A;) 



of a homogeneous system is proportional to the Fourier transform of the field correlation 
function g^^\x); a phase coherence extending over the whole sample therefore implies a 
macroscopic occupation of a single mode ^ . 

This finite-size induced Bose condensation occurs around the temperature 



^Apart from a numerical factor of the order of unity, the analogous condition for an interacting 
harmonically trapped cloud ic — Rtf {Rtf is the Thomas- Fermi radius of the cloud) leads to the 
same condensation temperature Tp^ = as predicted in from a different point of view. 
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at which the maximum number N^ax of atoms that can be stored in the excited modes 

1 mkBT 



ST^ 1 mkBTL"^ . 



k^O 

equals the number of atoms present in the sample; the longer L, the smaller the "tran- 
sition" temperature Tbec- In terms of the phase coherence length ^C1 the BEC condition 

> Nmax can be rewritten as L < 6£c- 
In fig.0 we have plotted the momentum distribution for a growing number of atoms across 
the transition value Nbec- for ^ Nbec the atoms that are further added to the system 
spread over the whole distribution. On the other hand, for N 3> Nbeci the population 
of the excited modes saturates to the value N^ax and the added atoms accumulate in the 
fundamental mode; this behavior is strictly analogous to what happens in three-dimensional 
systems. 

For the three upper curves for which N > 2Nbec, the Bogoliubov sampling technique of 



Sec.[illB| has revealed to be much more efficient than the Brownian motion one; for the lower 



ones, the Brownian motion technique of Sec.|TT^ was instead more performing. As a further 
check, we have successfully verified that the predictions of the two methods really agree with 
each other; on the scale of figJ3 the corresponding curves are nearly indistinguishable. 



2. Suppression of density fluctuations 

In addition to the features of the momentum distribution function and the related field 
correlation function g^^''{x) discussed in the previous section, the presence of a Bose- Einstein 
condensate can be detected also by looking at the density fiuctuations quantified by the value 
g^'^\0) of the second-order correlation function (0). For a non Bose-condensed gas, it is 
well-known that ^''•^■'(O) is close to 2 as for a chaotic beam of light; at zero temperature 
T = 0, all the A^ ^ 1 atoms of the system are in the fundamental mode, so that g^'^\0) = 
1 — ~ 1. As we immediately see in fig.H, the crossover between the two regimes occurs in 
the neighborhood of the Bose-Einstein condensation temperature Tbec- 

A simple argument allows to estimate the value of 5'*-^-*(0) at low temperatures, when the 
non-condensed fraction is small; the second order correlation function can in fact be written 



as 



\ k^k' k / 

V k^k' k k J 



The sum can be analytically computed since the system is assumed to be degenerate: in this 
case, most of the population is concentrated in the lowest modes for which j3ti?k'^ /2m <C 1, so that 
the exponential can be linearized e^^ ^ /^"^ — 1 ~ ph'^k'^ /2m. 
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Whenever the non- condensed fraction Nnc/N = 'Ylki'^k) / N is small enough, the probability 
for having a completely empty condensate results negligible; under this assumption, the 



non-symmetry-breeaking Bogoliubov approach of is exact for the non-interacting gas 
and the spectrum of excitations is composed of a set of harmonic oscillators, one for each 
7^ mode of the trap. At a given temperature T, each mode is thermally populated in an 
independent way with a mean occupation number 



{rik) 



and a variance 



By substituting ( |^ and (^) into (^Up, we get the simple result 



1 



1 + 2- 



N 



Ar2 

nc 



3E 



Ar2 



(51) 



(52) 



(53) 



which is plotted as a solid line in fig.^. The agreement with the prediction of the stochastic 
calculation is excellent especially at low temperatures, when the non-condensed fraction 
Nnc/N is substantially smaller than 1; the discrepancy becomes important only when the 
non-condensed fraction is as high as 1/2 and the probability of having an empty condensate 
mode is no longer negligible. 

It is worth stressing that the present suppression of density fluctuations is a consequence of 
the statistics only and therefore completely different from the one originating from interac- 
tions that we shall discuss in Sec. [IV CT] . 



3. Ground mode population statistics 

The probability of having n atoms in the condensate is given by the expectation value of the 
projector (5(o:n) on the subspace in which the ground k = mode contains exactly n atoms: 

Qo{n) = Tr[g(o:„)p] Tr[p]-' (54) 

In terms of our Hartree ansatz (^, this quantity can be written as 

A^l 

where 0° 2 components of the wave functions 0i^2 along the ground k = mode and 

(j)i2 components along the orthogonal subspace spanned by the remaining k ^ 

modes. In the non-interacting case in which (pi = (f)2 = <Pi this expression can be further 
simplified to 

A^l 
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At temperatures T much higher than the condensation temperature Tbeci the mean popu- 
lation of all the modes is only a small fraction of the total number of particles so that each 
of them sees all the others as a sort of particle reservoir and its occupation is accurately 
described within the grand-canonical ensemble by the thermal distribution function 

of average occupation number n. In this regime, the occupation probability Qo{n) of the 
ground mode /c = is a monotonically decreasing function and the most probable occupation 
number is n = (fig.^, solid line). 

As soon as T approaches the condensation temperature Tbec, the occupation statistics of 
the fundamental mode is drastically modified by the presence of a condensate (dotted and 
dot-dashed lines) and, in particular its slope at n = changes its sign around T = Tbec P^ - 
For T < Tbec the occupation function Qo{n) has a peaked shape around its mean value 
(dashed line). 

A similar transition is well known to occur for the photon distribution in a laser: below 
threshold the distribution has a thermal shape (^), while above threshold it tends to the 
Poisson distribution typical of a coherent state [pll]. 



C. Interacting gas 

In the case of the non-interacting system described in the previous section it has been possible 
to find the simple a priori distribution function P[0] (^) which optimizes the efficiency of 
the Monte Carlo calculation at a given /3. For an interacting system, such a optimum choice 
is no longer simple to implement; heuristically, the simulations of the interacting system 
described in the present section have been carried out using the same a priori P[4>] as used 
for non-interacting gases. 

For simplicity of the analysis, we shall limit ourselves to the case of a zero-range interaction 
potential 

Vint{x - x') = g 6{x - x') ; (58) 
the coupling constant g is related to the ID scattering length am by 

9 = (59) 

maiD 

On the discrete lattice we are actually using in the simulation, the delta interaction potential 
(^) translates into 

Vir,t{x ~ = ^ ^x,x', (60) 

Ax being the lattice spacing. 
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1. Classical field regime 



If the gas is well in the degenerate regime T <C Tdeg but not condensed, the main effect of 
atom-atom interactions is the suppression of density fluctuations and the appearance of a 
new length scale ^ = ^jWjlmng over which the density fluctuations are correlated ||2^, the 
so-called healing length. 

If interactions are extremely weak x = f?gfP"iv' jm <^ 1, i.e. ^ ^ £c, then the gas can be 
considered as an ideal gas, the density fluctuations are close to the non-interacting value 
(7*^^^(0) = 2 and the characteristic length scales for phase correlation and density fluctuations 
turn out being of the same order 4 (figia). The fact that c/(2)(o) in the fi gure is not exactly 
2 even in a non-interacting case is due to the fact that the system we are considering is 
a finite one; we have in fact checked that for increasing box size L, ^''■^-'(O) approaches its 
thermodynamical limit of 2 (fig.pD). 

For stronger interactions x > 1; density fluctuations are strongly suppressed by the interac- 
tions and their correlation is limited to the shorter length scale ^ < lc\ phase correlations 
are instead far less affected by the presence of the interactions than density correlations and, 
in particular, the field coherence length is increased by less than a factor 2. 

From a quantitative point of view, the agreement of the result of our stochastic simulations 
with the approximated predictions of classical field theory |2^ looks reasonable. We are 
in fact in a regime in which the classical field approximation can not be expected to be 
extremely accurate, since the classicity parameter ri = gn/ksT < 0.1 and the terms which 
were neglected are of order r]. 

If we neglect the stochastic noise in our calculations, we are led to an improved formulation 



of classical field theory: with respect to the usual one |24|, we are now including one more 



"force" term in the imaginary-time evolution of the wave function in addition to the classical 
Boltzmann weight e"^"^'"^'. The presence of this additional term makes the approach exact 
at least in the case of a non-interacting gas and, in particular, always immune from the 
so-called black-body catastrophe; the new deterministic term is in fact able to reproduce 
the correct Bose law even for the high energy modes, while the Boltzmann weight alone 
would predict a mean energy of ksT per mode. In practice, we have found that the result 
of the simulations is not appreciably modified if we neglect the stochastic term provided we 
are well within the classicity region 77 -C 1; in other terms, this means that the predictions 
of the improved classical field theory are accurate in the 77 -C 1 regime. 

A similar improved classical field theory can also be formulated within the coherent state 
approach if we neglect the non-positive diffusion term in the Fokker-Planck equation ||24|| ; 
also in this case, the approach gives exact predictions for the non-interacting gas. 
A characteristic feature of classical field theories based on the Glauber-P distribution such as 



the one in [24] is that the observables can be expressed as mean values of the corresponding 
classical field variables over a real and non-negative distribution function. Thanks to the 
Schwartz inequality, this implies that the second-order correlation g^'^\0) is always larger or 
equal to 1 and the field can never be antibunched. 

In the next section we shall push our simulations over the boundary 77 = 1 of the classic- 
ity region and we shall look for signatures of non-classical behavior such as antibunching 
(g^'^\0)). In previous stochastic field approaches such as the Glauber-P, such a parameter 
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regime could not be explored since the non-positive diffusion term could not be simulated 
by means of usual Monte Carlo techniques |pl|| . 



2. Beyond the classical field regime: spatial antibunching of atoms 



In the Sec.LV B we have seen that below the Bose-Einstein temperature Tbec a non- 



interacting gas is condensed in the zero-momentum mode. In the present section we discuss 
the effect of interactions on a condensed sample for t] > 1, i.e. outside the classicity region. 
Since the mean occupation number of the fundamental mode is large, Bogoliubov rejection 
sampling will be used for the stochastic simulations. 

In fig.^ we have plotted the value of ^''•^•'(O) as a function of temperature for different values 
of the interaction constant g; as T tends to zero, ^''•^•'(O) tends to a finite value which is 
equal to 1 — for the non- interacting gas and equal to a smaller value in the presence 
of repulsive interactions; the stronger the interactions, the lower g^'^\Q). In the other (b) 
panel, we show the behavior of g^'^\x) well below the condensation temperature: due to the 
repulsive interactions, a dip is present around x = 0. This effect is due to the stochastic 
term: the result of a calculation performed neglecting the noise term is shown as a long- 
dashed line; as expected g^'^\Q) tends to 1 — for T — > 0. The disagreement with the 
exact calculation was expected since the classical approximation can not be accurate in the 
present definitely quantum regime 77 ~ 5. 

We now wish to compare our results for g^'^\Q) to approximate analytical predictions. As 
a starting point we use the following result: given a Hamiltonian Ti which depends on 
a parameter g (in our case, the interaction coupling constant), it follows from elementary 
quantum statistical mechanics that the partial derivative of the free energy F = —ksT \ogZ 
[Z is the partition function Z = Tr[e~^'^^^'']) with respect to the interaction parameter g is 
equal to 

/^^\ /fill 

in the case of a local interaction potential (|58D , this brings to the remarkable expression 

which can be used to obtain an analytical prediction to be compared with the results of 
stochastic simulations. 

Furthermore, we are in a regime where the temperature is well below the transition tem- 
perature Tbec and the interactions are sufficiently weak n,^ 1, most of the particles are 
therefore in the condensate and the system is accurately described by a Bogoliubov approx- 
imation in which the Hamiltonian of the many-body system is approximated as a system of 
uncoupled harmonic oscillators corresponding to the different elementary excitations of the 
system Jll^. 



In this case, the free energy F of the system of uncoupled harmonic oscillators is immediately 
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calculated as 



F{g) = Eo{g) + keT log(l - e-"'") (63) 
and its derivative with respect to g has the following simple form 

dg dg f^^ dg e^'^ - 1 ^ ^ 

The index k ^ runs over the different modes whose energies are given by the well-known 
Bogoliubov expression 



2m \ 2m 

and the ground state energy -E'o(fl') is equal to |T 



^'^ = \h::r{^ + ^9n] (65) 



2L 

with Vk = Isinh and 



tanh 2ek = ■ (67) 

ek + gn 



Inserting the explicit expressions ( |65| , |66| , |67D into (|6^), numerically performing the sum over 
the different modes of a finite homogeneous box and finally inserting the result into (|62D 
we are led to the finite temperature Bogoliubov predictions for g^'^\0) which are plotted in 
fig.^a as dashed lines. 

These results can be compared to stochastic simulations: for T — > 0, Bogoliubov theory 
is most accurate for weak interactions, when the quantum depletion of the condensate is 
small. At finite temperatures (but much lower that the condensation temperature), thermal 
excitations are present on top of the condensate: since the free-energy ( |B3D only keeps track 
of the lowest order terms in the excitation density, the agreement of the approximated 
Bogoliubov predictions with stochastic simulations is best at low temperatures when the 
density of non-condensed atoms is small. 

In the non-interacting case, the stochastic simulations can be compared also with the theory 



of Sec. [lV B 2| ; as expected, that approach is more accurate that the Bogoliubov calculations 



of the present section (fig.^); in the Bogoliubov approach, in fact, the g = result is 
obtained as the limit for g ^ of a theory which is linearized in the excitation density and 
thus able to reproduce only the terms in (|53|) which are linear in the rifc's; the quadratic 
terms are instead not included at all. 



V. CONCLUSIONS 

In this paper, we have developed a new method for the exact calculation of the thermo- 
dynamical quantities of an interacting A^-boson system in the canonical ensemble. This 
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method is the imaginary-time version of the one discussed in our previous paper [jl^ for the 
real time evolution of the system. 

In order to have an efficient sampling of the observables, an importance sampling method 
has been used with an a priori distribution based on the non-interacting gas. We have 
applied the simulation scheme to a homogeneous one-dimensional interacting Bose gas and, 
whenever available, we have successfully compared the predictions of stochastic simulations 
with the ones of existing theories, such as the grand-canonical ensemble, the classical field 
approach [^] and Bogoliubov theory [12-15|. 

In particular, we have discussed the effect of Bose condensation in finite-size systems on the 
different observables, such as the momentum distribution, the density fluctuations and the 
occupation statistics of the ground mode. 

In the case of a non-condensed system, the main effect of interactions is to suppress density 
fluctuations without affecting in a dramatic way the phase coherence properties of the sam- 
ple. In the classical fleld regime ksT ^ /i, the effect of the noise term is negligible so that 
the system is already accurately described by the deterministic force term alone. For the 
opposite case of a condensed sample at very low temperatures ksT < /i, interesting features 
have been predicted such as an spatial antibunching of the atoms. 

In the future, we plan to combine the imaginary-time evolution discussed in the present paper 
with the real-time evolution discussed in in order to dispose of a tool for the numerical 
calculation of dynamical properties of flnite temperature Bose gases. Improvements of the 
statistical properties of the stochastic simulations will hopefully be investigated both in 
terms of better a priori distribution functions and in terms of the extension of the ansatz to 
include more sophisticated states such as Bogoliubov vacua. Finally, we plan to generalize 
the stochastic approaches in order to consider also multi-mode nonlinear optical systems in 



which driving and damping play an essential role [21 



The flrst part of the work on the real-time evolution has been done in collaboration with 
Jean Dalibard, to whom we are indebted for many stimulating discussions. Laboratoire 
Kastler Brossel is a Unite de Recherche de I'Ecole Normale Superieure et de I'Universite 
Paris 6, associee au CNRS. 
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FIG. 1. Homogeneous non-interacting gas in a large L ^ box with periodic boundary condi- 
tions. Mean occupation number of the different plane wave modes of the box in the non-degenerate 
(T = 8Tdeg) regime (a) and in the degenerate (T = 0.02Tdeg) regime (b); the disks are the result of 
stochastic calculation with = 2048 realizations, the solid line is the prediction of grand-canonical 
ensemble in the thermodynamical limit. = 6 atoms (a) in a box of length L = 6 = 96 atoms 
(b) in a box of length L = 961 {I is an arbitrary length unit). 




FIG. 2. Homogeneous non-interacting gas in a large L ^ ic box. First (a) and second (b) order 
correlation functions below degeneracy for T = 0.02Tcieg (dot-dashed) and T = 0.08 T^gg (solid). 

= 96 atoms in a box of length L = 96 / (/is an arbitrary length unit). Stochastic simulation with 
J\fr = 400 realizations. The deviation of (7^^^(0) from the grand-canonical prediction g^'^\0) = 2 is 
a finite size effect (see Sec. IV B 2.) 



22 




FIG. 3. Degenerate non- interacting gas in a finite box: mean occupation number of tlie different 
momentum modes for growing number of particles across the Bose-Einstein condensation threshold 
N/Nbec = 0.25 (circles), N/Nbec = 0.5 (squares), N/Nbec = 1 (diamonds), 2 (triangles), 3 
(crosses) and 4 (stars). Stochastic calculation with Mr = 2048 realizations. L = 48/, /? = 16 {I 
is an arbitrary length unit). 




FIG. 4. Degenerate non-interacting gas in a finite box: second-order correlation function across 
the Bose-Einstein condensation threshold at T ~ Tbec- Solid line: analytical prediction (|53|). 
Disks: stochastic calculation with Mr = 256 realizations. = 192 atoms in a L = 48/ box (/ is an 
arbitrary length unit). 
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FIG. 5. Non-interacting Bose gas: probability distribution of the ground mode occupation 
number for different temperatures across the Bose-Einstein transition: T/Tbec = 30 (sohd), 3 
(dotted), 1 (dot-dashed), 0.5 (dashed). Stochastic calculation with AT^ = 1024 realizations. N = 32 
atoms in a L = 24Z box {I is an arbitrary length unit). 
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FIG. 6. Left panel: first g^^\x) and second-order g^'^\x) correlation functions for different 
values of the interaction strength g/{f? /ml) = (solid), 0.005 (dotted), 0.01 (long-dashed) and 0.02 
(dot-dashed) (% = h'^gP'^n^ /m = 0, 1, 2, 4), all within the classical field regime [gn/kBT < 0.1 ^ 1). 
Right panel: comparison of the approximate classical field predictions in the thcrmodynamical 
limit with stochastic simulations for finite boxes of growing size L/ic = 9 (circles), 18 (squares). 
Stochastic simulation with J\fr = 2048 realizations, /3 = 0.8^. N = 324, L = 48/ (circles); 
N = 648, L = 961 (squares) {I is an arbitrary length unit). 
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FIG. 7. Left panel: second order correlation g^'^\0) as a function of temperature for different 
values of the interaction strength; from below = 0.1,0.05,0. Disks: stochastic calculation 
with A^r- = 2048 realizations. Dashed line: Bogoliubov prediction ( |62| , |6^ ) . Dotted line: analytical 
prediction ( [5^ ) for the non-interacting gas (nearly superimposed to the disks). Dot-dashed line: 
stochastic calculation neglecting the noise term for = 0.05; notice the low temperature limit 
equal to 1 — 1/N. Right panel: second order correlation function g^'^\x) for ^g = 0.1 (solid) at 
Tbec/T = 21 obtained by means of a stochastic calculation with Mr = 2048 realizations. = 42 
atoms in a -L = 6 / box (I is an arbitrary length unit). 
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